Layering, freezing and re-entrant melting of hard spheres in soft confinement 
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Confinement can have a dramatic effect on the behavior of all sorts of particulate systems and 
it therefore is an important phenomenon in many different areas of physics and technology. Here, 
we investigate the role played by the softness of the confining potential. Using grand canonical 
Monte Carlo simulations, we determine the phase diagram of three-dimensional hard spheres that 
in one dimension are constrained to a plane by a harmonic potential. The phase behavior depends 
strongly on the density and on the stiffness of the harmonic confinement. Whilst we find the familiar 
sequence of confined hexagonal and square-symmetric packings, we do not observe any of the usual 
intervening ordered phases. Instead, the system phase separates under strong confinement, or forms 
a layered re-entrant liquid phase under weaker confinement. It is plausible that this behavior is 
due to the larger positional freedom in a soft confining potential and to the contribution that the 
confinement energy makes to the total free energy. The fact that specific structures can be induced 
or suppressed by simply changing the confinement conditions (e.g. in a dielectrophoretic trap) is 
important for applications that involve self- assembled structures of colloidal particles. 

PACS numbers: 64.75.-g, 68.65.Ac, 64.70.-p, 82.70.-y 



INTRODUCTION 



The behavior of particles in confined geometries is im- 
portant in many different areas of physics and technology. 
This includes the physics of ions in electromagnetic traps 
[l| , of dusty plasmas confined by external fields 0, Q , of 
classical electrons in quantum wells [4J and of colloidal 
suspensions in narrow slits [5[, as well as application- 
oriented topics in nanotechnology, (bio) lubrication and 
the self-assembly of microstructured materials (e.g. Ref. 
(6l-[lQj). It is well-known that confinement effects can 
dramatically change the behavior of such systems, both 
quantitatively and qualitatively. For instance, when a 
suspension of colloidal hard spheres is confined in a 
wedge-shaped geometry, one observes a rich cascade of 
different equilibrium crystal structures as the wall spac- 
ing increases pj14T3| . Such behavior contrasts sharply 
with the bulk phase diagram, which consists of a sin- 
gle, density-dependent liquid to face-centered-cubic crys- 
tal transition. The past few decades have seen many 
studies of the behavior of hard spheres between impene- 
trable walls (e.g. Ref. [l3l4l6j), of charged particles un- 
der strong confinement (a model for trapped Coulombic 
or Yukawa particles, e.g. Ref. B TTtL fig ) , and of con- 
fined dipolar colloids (e.g. Ref. |19l-l2l|). However, to the 
best of our knowledge, there have been no studies that 
systematically investigate the role played by the softness 



of the confining potential. This is an important issue to 
address, as more and more systems become available that 
involve some form of soft confinement. In the areas of 
colloid science and nanotechnology one can for instance 
think of suspensions confined in dielectrophoretic (22I - 
[2^ or laser-optical fields 0, |2fj, nanometric objects in 
charged slits j gj, p articles interacting with soft polymer 
substrates [23, [28| , or particles trapped at liquid-liquid 
or liquid-gas interfaces [29|-[33|. Using Monte Carlo simu- 
lations, we here study three-dimensional systems of hard 
spheres that in one dimension are constrained to a plane 
by a harmonic potential. Unlike hard boundaries, this 
soft potential well does not prescribe a particular con- 
finement 'width', can be continuously tuned from very 
strong to very weak confinement, and makes an energy 
contribution to the total free energy which depends on 
the exact particle positions. We highlight the unique 
properties imparted by such soft confinement by com- 
paring the observed phase behavior with that of hard 
spheres between two hard walls [14[, as well as with the 
behavior of a more complex system of highly charged par- 
ticles and their counterions between neutral walls [34[ . In 
the latter system, the particles experience a combination 
of an effective harmonic potential due to the counterions 
and long-ranged repulsive Coulomb interactions. 



II. MODEL 



We performed grand canonical Monte Carlo simula- 
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each other unless their cores overlap and for any pair of 
particles the interaction potential is given by: 

^sphere— sphere / \ J 0, T ^ (J 

k^f ~ \ oo, r < a 

where ks is the Boltzmann constant, T is the absolute 
temperature, r is the center-to-center distance of the par- 
ticles and <j is the particle diameter, which we took as 
the unit of length in our simulations. We used periodic 
boundary conditions in the xy plane (box size L x x L y 
= 20 x 20, unless stated otherwise) and a soft confining 
harmonic potential centered around z = 0, which acts 
on each particle individually and whose softness was set 
through the spring constant k (in units of <r -2 ): 



We performed simulations for different chemical poten- 
tials of the reservoir, starting with an empty box. Using 
cell lists and an early rejection scheme [35], we typically 
performed 6 x 10 10 Monte Carlo moves in which we at- 
tempted to insert, delete or displace (d max = 0.05) a 
randomly chosen particle, where the fraction of insertion 
and deletion moves was fixed at 0.2. We also performed 
simulations with an additional move in which the L x /L y 
aspect ratio of the box was allowed to change while keep- 
ing the area A = L x L y constant. In principle, this should 
prevent the box shape from dictating the structure of the 
particle packing, but in all cases the resulting structures 
were identical to those observed in a square box. 

After thorough equilibration, we determined the num- 
ber of particles N in the simulation box and from this 
the density, p = N/A, and characterized the structure 
of the typically layered particle packings by calculat- 
ing the four- fold (q^) and six- fold (q^) symmetric two- 
dimensional bond order parameters in each of the lay- 
ers [36[. The bond-order parameters consider all of the 
nearest neighbors of a given particle that lie approxi- 
mately in the same z plane. Based on geometric ar- 
guments and the observed particle distributions along 
the z axis, the nearest neighbors were here defined as 
those particles residing within a center-to-center dis- 
tance tnn ^ 1-3 of the particle of interest and with a 
height difference £)znn ^ 0.3. This definition excludes 
second- nearest neighbors, which in a close-packed square- 
symmetric layer reside at r = V% and is also found to 
reliably discriminate between thermally broadened and 
adjacent layers along z. The ratio between the two 
bond order parameters allowed us to distinguish between 
the three main types of structures found in our simula- 
tions: disordered liquid for 1/3 ^ q±/ 1 q§ ^ 3, hexagonally 
(A) packed for q^/q^ < 1/3, and square (□) packed for 
#4/^6 > 3. The crossover values were determined from 
plots of q±/q$ across the entire density range studied in 
the simulations, as well as the density probability distri- 
butions, which revealed the first-order phase transitions. 
We find that varying tnn and £znn by ±0.1 does not 
affect the resulting phase diagram. 



III. RESULTS 

The phase diagram in Fig. [I] and the simulation snap- 
shots in Fig. [2] provide an overview of the observed parti- 
cle packings as a function of the density and spring con- 
stant. At high spring constants (k > 150), i.e. relatively 
strong confinement, and low densities the particles form 
a layer of disordered liquid (L) at the minimum of the 
confining harmonic potential. At higher densities, the liq- 
uid freezes into one or more crystalline layers which have 
a hexagonal (A) or square (□) symmetry (Fig. [2]), fol- 
lowing an alternating sequence as the density increases: 
1A 2D 2A -> 3D 3A . . . (the integers indi- 
cate the number of particle layers that is formed) . To ra- 
tionalize these observations, we consider the limit k — » 00 
(or T — »• 0) in which the entropic contribution to the free 
energy can be neglected and the energy contribution of 
the confining potential dominates: F « E = fj] 2 ^ 
where the sum runs over all the particles in the system. 
Simple geometry then gives the energies of the perfect 
hexagonal and square packed structures, which scale lin- 
early with the density because of the quadratic form of 
the confining potential, Fig. [3l The maximum density 
of each of the phases corresponds to close packing, e.g. 
Pmax = 2 for 2D and p max = 4/V^ for 2 A. Starting at 
low density, the L — >> 1A transition is given by the two- 
dimensional hard disk freezing transition [37], [HI and the 
1A phase (strongly constrained to z = and with neg- 
ligible energy) is stable until close packing at p = 2/^/3- 
Instead of continuously transforming into a 3 A structure, 
the system then phase separates into 1A and 2D. This 
phase separation can be understood in terms of a free 
energy minimization criterion, equivalent to the double- 
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FIG. 1: (Color Online) Phase diagram in the spring constant - 
density representation. Every point represents the result of a 
single simulation run. Black crosses ( x ) : liquid. Red triangles 
(A): hexagonally packed. Blue squares (□): square packed 
(structures shown in Fig. [2]). Approximate phase boundaries 
are indicated with solid lines. Dashed lines indicate the stable 
phases in the limit k — >■ 00 (from Fig. [3]). 
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FIG. 2: (Color Online) Simulation snapshots at different den- 
sities (for k = 100) and 'ideal' schematic representations 
(lower right) of the spatial arrangement of the particles in 
the square (□) and hexagonally (A) packed structures. 

tangent construction. A similar argument applies to the 
subsequent transitions 2D — >• 2 A, 2 A — » 3D, and so on. 
Around 7 layers the square symmetric structures eventu- 
ally disappear, in favor of the denser hexagonal packings, 
which have a face-centered cubic (fee), hexagonal close 
packed (hep) or a random hexagonal close packed (rhep) 
structure, similar to the crystalline bulk phases of hard 
spheres. 

The sequence of alternating hexagonal and square 
packings under strong harmonic confinement corresponds 
to the simple Pieranski picture of hard spheres con- 
fined between two hard walls [11]. However, we do not 
see any of the buckled, rhombic and prism phases that 
are found to interpolate between these packings under 
hard- confinement conditions, when entropy solely deter- 
mines the phase behavior [12l-ll4j. By contrast, at high 
spring constants the behavior of the harmonically con- 
fined system is energy-dominated and the usual interven- 
ing phases are found to be unstable, because the second 
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FIG. 3: (Color Online) Normalized energies per unit area of 
the hexagonal (red lines, A) and square (blue lines, □) packed 
structures in the limit k — >> oo. 



derivative of the (free) energy with respect to the den- 
sity is negative. Furthermore, the intervening phases also 
appear to be unstable at finite spring constants, as we 
always observed a spontaneous melting (low spring con- 
stants) or phase separation (higher spring constants) of 
the system when it was initially prepared in one of these 
phases. In the absence of any stable intervening phases at 
higher spring constants, the coexistence regions between 
the stable hexagonal and square packed phases - which 
can only exist with low free energy at certain densities 
due to the integer number of layers - are wider than under 
hard-confinement conditions (note, by the way, that in 
order to approach the hard-confinement limit one would 
need to increase the exponent of the confining potential, 
rather than the pre- factor). We point out that the alter- 
nating hexagonal and square-symmetric packings appear 
to be a common property of different types of confined 
repulsive particle systems [5|, LLJ, [34], [39| , while the char- 
acter of the intervening phases seems to depend more 
strongly on the exact details of the particle-particle in- 
teractions and the confining potential. For example, at 
high spring constants we do not observe any new phases 
between 1A and 2D in the harmonic potential, while the 
same hard spheres in a hard slit would form an interme- 
diate buckled bilayer structure (223) [l2l-[l4]], and while 
charged particles in an effective harmonic potential are 
expected to form the sequence 1A — >• 3 A — >> 2D in the 
limit T — > [34|]. In the latter system, which considered 
point-like particles, the long-ranged repulsive Coulomb 
interactions between the particles compete with the at- 
traction to the minimum of the confining potential. 

Remarkably, as the harmonic confinement becomes 
softer (lower spring constant), we do find stable inter- 
vening phases, which, however, are not ordered. Instead 
of the more commonly observed solid-to-solid transfor- 
mations, the system undergoes a couple of re-entrant 
melting transitions that give rise to intervening liquid 
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FIG. 4: (Color Online) Density-dependent re-entrant melting 
between the 1 A and 2D phases, for k = 100 and different sizes 
of the simulation box (L = L x — L y ). Red stars: L — 50. 
Green diamonds: L — 20. Blue circles: L — 10. The insets 
give the probability distribution of observing a given density 
at coexistence chemical potential on the same density scale as 
the main plot. 



phases, with triple points around {k = 120, p = 1.2} 
and {k = 45, p = 2.25}. Thus, around k « 100 we 
observe freezing into a stable 1A phase, which then re- 
melts into a disordered liquid before freezing again into 
the 2D phase (Fig. [2j), while for higher densities we find 
the same alternating sequence of hexagonal and square 
packings as before. We further see that the softer the 
confinement, the higher the density at which the initial 
liquid phase still persists and at sufficiently low spring 
constants (k < 80) the first stable ordered structures ac- 
tually consist of more than one layer. For example, at 
k « 40 the first ordered phase is 2 A, followed by re- 
entrant melting, the 3D structure, and then the other 
multi-layered ordered phases, whereas for k < 35 we only 
observe the latter, without any intermediate melting. 

To make sure that the re-entrant liquid phases did not 
represent finite-size artifacts, we performed simulations 
with three different box sizes (L x = L y = 10, 20 or 
50) for k = 100, Fig. 2J It can be seen that the re- 
sults for the different system sizes are essentially iden- 
tical and that there is a clear transition from the 1A 
phase to a re-entrant liquid and then to the 2D phase, as 
reflected by the q^/q^ bond order parameter ratio. We 
also determined the density probability distributions at 
coexistence [5o|, |4l| using histogram-reweighing [42| and 
expanded ensemble [43j simulations (insets in Fig. |4j). 
For both the 1A — >> L and the L — >• 2D transitions we 
find bimodal probability distributions that are charac- 
teristic of a first-order phase transition. In addition, Fig. 
[5] shows the spatial distribution of the particles in the 
soft confining potential for different densities outside the 
coexistence region. As expected, in the 1A phase the 
particle distribution has a single peak centered at the 
energy minimum z = and in the 2D phase there are 
two peaks which are symmetrically located with respect 




z 



FIG. 5: (Color Online) Particle distribution along the axis 
of confinement (z), for k = 100 and different densities. Red 
triangles: 1A, p — 1.05. Green crosses: liquid, p — 1.23. 
Yellow diamonds: liquid, p = 1.43. Blue squares: 2D, p = 
1.58. 

to z = 0. The intervening re-entrant liquid starts out 
with a flat and broad distribution of the particles at low 
densities, which then continuously transforms into a dou- 
bly peaked profile as the density increases. Interestingly, 
this occurs without further ordering of the particles in 
the xy plane, so that the result is a layered liquid. Oth- 
ers have shown that if the same hard spheres are confined 
between impenetrable hard walls, a highly ordered buck- 
led structure forms when the wall separation is larger 
than required for the 1A phase, but too small for the 
2D phase [l2l-[l4j]. This buckled 2B phase optimizes the 
particle packing in the available space between the two 
walls by splitting the 1A structure into rows of particles 
that alternate in height, effectively forming a crystal of 
two interpenetrating rectangular layers, which continu- 
ously transforms into the 2D structure as the wall sepa- 
ration increases. In the harmonically constrained system, 
the re-entrant layered liquid fulfills a similar interpolat- 
ing role between the one and two-layer ordered phases, 
but the additional energy considerations and greater po- 
sitional freedom under soft confinement, together with 
the absence of long-ranged particle interactions favor dis- 
ordered over ordered phases. 



IV. CONCLUSIONS 

We have identified the stable phases of hard spheres 
under harmonic confinement, as a function of the den- 
sity and the softness of the confining potential. We find 
the well-known 'base' sequence of alternating hexagonal 
and square-symmetric packings (1A — >> 2D — >> 2A — >• 
3D — » 3 A — )> ...), which thus appears to be quite in- 
sensitive to the details of the (confining) interactions. 
When we look at the presence of stable intervening phases 
the soft-confined hard-sphere system is clearly differ- 
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ent from hard- confined systems, though. Instead of the 
usual highly ordered interpolating particle packings, we 
observe phase separation between the hexagonally and 
square packed structures under strong confinement, and 
disordered re-entrant liquid (L) phases under weaker con- 
finement. For 1A — » L and L — » 2D, we have shown 
that the re-entrant melting/freezing transition has a first- 
order character and that the liquid develops a layered 
structure, without further ordering of the particles. We 
argue that the re-entrant and phase separating behav- 
iors are due to the fact that the penetrable soft har- 
monic potential on the one hand offers more positional 
freedom than hard- wall confinement does, while, on the 
other hand, making an important energy contribution 
to the overall free energy of the system which depends 
on the exact particle positions (in addition to the usual 
entropic considerations for hard spheres). The net re- 
sult is a smaller diversity of crystalline packings, as com- 
pared to hard spheres between two hard walls. We expect 
that other soft confining potentials may well have a sim- 
ilar effect, although the details of the phase diagram will 
likely be different. The fact that certain ordered struc- 
tures can be reliably obtained, while many other struc- 
tures can be induced or suppressed on demand through 



a variation of the confinement conditions is important 
for applications that involve self-assembled structures of 
colloidal particles. Interestingly, it should be possible to 
realize harmonic confining potentials with dynamically 
tunable softness experimentally, for instance in a dielec- 
trophoretic trap. 
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